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Abstract — We use the quasi-static particle-in-cell code Quick- 
PIC to perform full-scale, one-to-one LWFA numerical experi- 
ments, with parameters that closely follow current experimental 
conditions. The propagation of state-of-the-art laser pulses in 
both preformed and uniform plasma channels is examined. We 
show that the presence of the channel is important whenever 
the laser self-modulations do not dominate the propagation. 
We examine the acceleration of an externally injected electron 
beam in the wake generated by ~ 10 J laser pulses, showing 
that by using ten-centimeter-scale plasma channels it is possible 
to accelerate electrons to more than 4 GeV. A comparison 
between QuickPIC and 2D OSIRIS is provided. Good qualitative 
agreement between the two codes is found, but the 2D full 
PIC simulations fail to predict the correct laser and wakefield 
amplitudes. 

Index Terms — Simulation, Lasers, Plasma waves, Electron 
accelerators 



I. Introduction 

The advent of high intensity (I > 10 18 W/cm 2 ), tightly 
focused (Wo ~ 10/irn) and ultra-short (tq > 30fs) laser 
pulses lead to the production of quasi-mono-energetic electron 
bunches via Laser Wake Field Acceleration (LWFA) [l]-[5]. 
These lasers are used in state-of-the-art experiments, leading 
to the self-injection of plasma electrons once the 'bubble' or 
'blowout' regime [6]-[8] is reached. Ideally in this regime, 
the transverse laser ponderomotive force is strong enough 
to expel nearly all electrons away from the laser axis. As 
electrons return to the axis they cross and form a thin electron 
layer, leaving behind an ion column. The spherical shape of 
the wake has ideal properties for electron acceleration, with 
linear accelerating and focusing forces [7], which provides 
high-energy and high-quality electron acceleration. 

The physics associated with the blowout regime is complex 
and can only be fully explored using numerical simulations, 
which constitute an insightful tool to examine the mechanisms 
associated with electron acceleration in this regime. The 
increased interaction length required for higher energy gains 
poses, however, a natural difficulty to full PIC simulations, 
which require more than 10 5 CPU hours for the model- 
ing of centimeter-scale plasmas; therefore systematic studies, 
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modeling and planning of new LWFA experiments with full 
PIC simulations, become increasingly difficult. One way to 
overcome this difficulty is the use of reduced PIC codes, 
such as QuickPIC [9] which allow for computational time 
savings of more than three orders of magnitude in comparison 
to full PIC codes. In this paper, we will then use QuickPIC 
to perform one-to-one, full scale LWFA simulations, both in 
uniform plasmas and in plasma channels, with and without 
externally injected electron beams. 

This paper is structured as follows. In section III] we briefly 
present QuickPIC physical model. In Section |Hl| we examine 
~ 2 J laser evolution, and the corresponding plasma response, 
using physical parameters that closely follow Refs. [4], [5]. 
In addition we present simulation results with high energy 
laser pulses, over longer propagation distances, towards stable 
multi-GeV electron acceleration. In Section |TVl a critical 

and we 



discussion of the results of Section III 



is given, 



compare both the laser evolution and driven wakefield between 
QuickPIC with full 2D OSIRIS [10] simulations. Finally, we 
state the conclusions. 

II. QuickPIC physical model 

For the sake of completeness we outline here the key 
assumptions of the QuickPIC physical model [9]. 

QuickPIC works under the quasi-static approximation 
(QSA) [11] which is valid whenever the typical evolution 
time/length for the driver is much larger than the typical 
time/length for the plasma response. The typical time for the 
evolution of a laser driver in the QSA is generally determined 
by the Rayleigh length, Z r , while the plasma response is 
determined by the plasma wavelength, \ p . Thus, the QSA 
is always satisfied for the typical experimental parameters 
required for LWFA. The QSA fails, however, to model trapped 
particles, precluding the physics that is present in the self- 
injection. Consequently, the dynamics of electrons in the 
blowout regime can only be examined through the evolution 
of an externally injected electron beam. 

QuickPIC solves the Maxwell's equations, written in the 
Lorentz gauge under the QSA in the co-moving frame 
(x,y,£ = z — ct, s = z) given by: 



= 4^o, 
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where <f> and A represent the plasma scalar and vector poten- 
tials, p and j represent the charge and current distributions and 
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= d 2 /d 2 x + d 2 /d 2 y. The charge distribution is initialized 

as: 

(An r 2 \ 

falling linearly for r > r c . The background plasma density is 
denoted by n , the depth of the channel by An, the distance to 
the axis by r = \J x 2 + y 2 and r c is the radius of the channel. 

Eqs. ([T| and |2]) reveal that the plasma response is fully 
determined by the transverse radial charge and current dis- 
tributions. Therefore, the simulation domain is divided into 
transverse slices, perpendicular to the laser propagation di- 
rection, where the plasma response is determined for a given 
laser profile. After solving the plasma response for the whole 
simulation box, the laser is advanced according to the pon- 
deromotive guiding center approximation: 

where A lascl is the laser vector potential, defined according 
to: 

a = a (I0r 3 - 15r 4 + 6r 5 ) exp [-r 2 /W 2 ] , (5) 

where a = qA lascT /(m e c 2 ), r = s/tq, tq is the laser duration 
at full width half maximum on the fields, and Wo is the spot 
of the laser. 

Even though Eq. Q evolves the envelope of the laser, 
the longitudinal grid size A£ should be chosen as small 
as possible, such that it is guaranteed that the frequency 
shifts associated with the envelope modulations are captured 
correctly. Since the discretization of Eq. |4]) implies a stability 
condition for the longitudinal grid cell size, given by A£ c > 
1/fco, A£ should be chosen as close as possible to A£ c . On 
the other hand, since Eq. Q must also resolve correctly the 
Rayleigh length, the spacing between successive s, As, must 
be chosen appropriately. 

QuickPIC has been benchmarked with full 3D OSIRIS 
simulations for a number of situations, both for the plasma 
wakefield accelerator [9] and LWFA [12]. Excellent agreement 
between the two codes was found. 

III. One-to-one modeling of state-of-the-art 

EXPERIMENTAL CONDITIONS 
A. Modeling of recent LWFA experiments 

We now examine the propagation of ~ 2 J laser pulse using 
state-of-the-art experimental parameters [4], corresponding to 
the evolution of a 12 TW and 25 TW laser pulses, in both 
uniform and preformed parabolic plasma channels. We use 
a Q = 0.8 and a = 1.1 for the 12 TW and 25 TW laser 
pulses respectively. In both simulations, the laser interacts with 
a plasma with uq = 3.4 x 10 18 cm~ 3 , and with An = 0.9 x 
10 17 cm~ 3 and r c = 110 zxm, if the channel is present. The 
laser pulse duration is t = 73 fs, focused to Wq = 25 /im. 
The simulation box is (360 /im) 2 wide and 60 [im long, 
divided into 512 2 x 256 cells for the transverse and longitudinal 
directions, respectively, with 4 particles per cell. 

We first examine the laser evolution of the "lower" power 
pulse in the plasma channel. Initially, the laser slightly disturbs 



the plasma, creating a perfectly linear wake with an amplitude 
well below the threshold for wave-breaking and self-injection 
(Fig. [T]-a). After entering the plasma, the pulse is strongly 
self-focused leading to a fast compression of the spot-size to 
roughly 15 /tm after only 2.5 mm of propagation. Simulta- 
neously, the peak vector potential of the laser increased to 
ao ~ 2 (Fig. fTJ-b) leading to an increase of the peak power 
of the pulse by nearly a factor of 2. In addition to self- 
focusing, this also resulted from longitudinal self-modulation 
effects [13]. As the laser is modulated and the transverse and 
longitudinal ponderomotive forces become stronger, electrons 
are progressively expelled from the axis, leading to nearly full 
electron cavitation at s sa 2 mm. By then, the accelerating 
gradient attains its maximum value, as shown in Fig. [TJ-a, 
suggesting that if self-injection occurs, it would be held most 
likely near that region. 

For s > 2 mm, the laser remains guided with spot 
oscillations of roughly 20% with a wavelength As = 2 mm. 
In order to understand if the guiding of the laser is only due 
to self-focusing, or due to the presence of the channel, we 
have performed a similar simulation in an uniform plasma. The 
simulation reveals that for s > 2 mm the laser continuously 
diffracts, thus making the channel of fundamental importance 
for the guiding of the laser. 

The main features associated with the 12 TW laser pulse 
propagation in the channel are also present in the 25 TW 
laser simulation, also in a channel, as seen in Fig. [2] De- 
spite being more intense, the initial wake is still below the 
threshold for wave-breaking and self-injection. However, as 
the laser propagates, the wake becomes stronger since both 
the laser spot is compressed due to self-focusing, and the 
peak vector potential is increased due to a combination of 
self-focusing, self-steepening and self-compression (Fig. |2ja). 
For s — 1.5 mm, the spot is compressed to its minimum 
value of 14 fim, and nearly all electrons are expelled from the 
laser axis. This leads to the excitation of a strongly non-linear 
wake (Fig. |5j-a), where the occurrence of wave-breaking and 
self-injection becomes likely. 

A distinct behavior, in comparison to the "lower" power 
laser scenario, occurs for s > 1.5 mm, where the peak vector 
potential keeps rising (Fig. |2]-b). This behavior is associated 
with the higher initial o,q of the 25 TW laser, which further 
enhanced longitudinal self-modulations. In order to understand 
the role of the channel in this mechanism, an identical simula- 
tion in an uniform plasma was performed, which revealed that 
the laser pulse remained guided even for s > 1.5 mm (Fig. (|2]- 
c)). Thus, both longitudinal self-modulations, in particular self- 
compression which increases the peak laser vector potential, 
and transverse modulations which are further enhanced by the 
latter, are sufficient to reach stable laser propagation. We note, 
however, that for propagation distances larger than those we 
have modeled (s max ~ 1 cm), the laser may not be self-guided 
as the transverse profile of the laser may develop higher order 
gaussian beams. Furthermore, it is important to observe that 
in the 25 TW scenario the self-guiding is directly associated 
with the structure of the bubble [7], [14]. 
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B. Multi-GeV electron acceleration in the blowout regime 

Until now, the LWFA experimental results have been ob- 
tained in the self-injection regime. However, the external 
injection of an electron beam to be accelerated by a laser 
driven plasma wave, constitutes a rather promising path to pro- 
duce quasi-monoenergetic, multi-GeV electron beams, without 
resorting to self-injection, because both the quality of the 
injected beam and the optimal accelerating properties of the 
blowout region can be carefully controlled. In this Section we 
then examine the dynamics of an externally injected beam, 
loaded in the bubble, driven by a 10 J and 15 J laser pulse. 

The externally injected electron beam is defined by a bi- 
gaussian distribution: 



n b = n h0 exp(-r 2 /2cr 2 )exp(-£ 2 /2crf), 



(6) 



where n^o is the peak electron beam density, ay and are 
the transverse width and duration of the beam. We note that 
QuickPIC treats self-consistently the fields associated with the 
external electron beam. 

We start by examining the propagation of a 10 J laser pulse, 
following [15], in an uniform plasma with n = 10 18 cm -3 . 
The laser is focused to Wo = 34 /im, has a duration of 
To = 30 is corresponding to a peak vector potential ao — 3. 
The beam is placed close to the rear part of the first plasma 
wave, in the region of the highest accelerating fields. The total 
charge of the beam is 11.2 pC, with ay = = 1.5 /im, 
roughly corresponding to the dimensions of a self-injected 
bunch in similar conditions. The beam charge was chosen to be 
far from the beam loading limit. The beam was injected with 
7 = 40, higher than the relativistic factor associated with the 
plasma wave phase velocity, which guarantees that the beam 
is trapped from the beginning of the simulation. The laser 
enters a simulation box with (330 /im) 2 wide and 60 fjm 
long, divided into 128 3 cells, with 4 particles per cell. 



In contrast with the simulations presented in Sec. III-A 



plasma electrons are almost fully blown-out from the laser 
axis right at the entrance of the plasma, creating a plasma 
channel which self-guides the laser with spot variations on the 
order of 30%. These transverse oscillations, however, do not 
interfere with the acceleration process, as beam electrons are 
accelerated with a constant acceleration gradient of roughly 
« 0.85 GeV/cm. 

An important property of the electron beam is the energy 
spread, which should be as small as possible. Since the electric 
field is stronger near the back of the bubble, beam electrons 
which are closer to that location accelerate faster. Naturally, 
this leads to an energy chirp and to the variation of the energy 
spread of the beam, which grows from AE lmtml / E < 1% to 
AE mm /B « 20%, at FWHM, until s < 1 cm. 

For s > 1 cm the opposite effect occurs. As the laser 
becomes pump depleted, the plasma wavelength decreases in 
such way that the front of the beam accelerates faster than 
the back. This effect progressively balances the beam energy 
chirp acquired until s < 1 cm, while it reduces its energy 
spread (phase-space rotation) thus leading to the formation of 
a quasi-monoenergetic bunch. After s = 2.3 cm, the beam 
then becomes nearly mono-energetic, with AE Rnal /E w 4% 



and with E as 2.2 GeV, as shown in Figure [3] We note 
that this phase-space rotation mechanism is different from that 
of Ref. [16]. In that case, the generation of the quasi-mono- 
energetic electrons was due to the dephasing of the beam in the 
plasma wave while in the present situation, the phase rotation 
of the beam is due to the pump depletion of the laser. 

It is possible to obtain an extremely stable acceleration 
regime, with negligible laser and wake oscillations, by using 
a preformed plasma channel and by carefully choosing the re- 
maining laser parameters. We have examined the acceleration 
of an externally injected beam in the wake of a 15 J laser pulse, 
with the laser and plasma parameters chosen according to the 
scaling laws to the blowout regime [7], [14]. The most stable 
laser propagation is obtained if the laser spot size is matched 
to the blowout radius of the electrons, rh w 2yfaQc/uj p , which 
leads to Wq ~ ~ 37/im, for ao = 2. In order to optimize the 
energy transfer from the laser to the beam electrons, the pump 
depletion length is matched to the dephasing length, which 
leads to To = 87 fs. The above considerations also imply that 
the background density is given by no = 1.6 x 10 17 cm~ 3 . 
We use a channel depth An = 6 x 10 16 cm~ 3 , which is 
3/4 of the linear guiding theory matching condition for a 
parabolic channel. We use an externally injected electron beam 
with 11.2 pC, chosen to be far from the beam loading limit, 
placed in the region of maximum accelerating fields, with 
ay = <T£ = 3 /im, roughly corresponding to the dimensions 
of a self-injected bunch in similar conditions. The beam is 
injected with 7 = 200, which is higher than the relativistic 
factor associated with the wave phase velocity, as determined 
from [14], [17] guaranteeing that particles are immediately 
trapped at the beginning of the propagation. The laser enters 
the simulation box with (612 urn) 2 wide, 140 /im long, 
divided into 128 2 x 256 cells for the transverse and longitudinal 
directions respectively, with 4 particles per cell. 

These parameters provide a nearly non-evolving acceler- 
ation structure, where the laser spot size oscillates within 
only 3% of the initial value. As the laser enters the plasma, 
it creates a mildly nonlinear wakefield, where the beam 
electrons dephase in the presence of an accelerating gradient 
of 200 MeV/cm, for the peak of the energy distribution. 
This is in good agreement with the theoretical predictions [7], 
[14] which predict an accelerating gradient of 260 MeV/cm. 
At the end of the simulation, after 20 cm of propagation, 
beam electrons gained more than 4 GeV. The theoretical 
maximum energy gain is sa 7 GeV after 26 cm, but for 
such long distances the laser envelope modulations are too 
strong, and the envelope approximation, critical for the laser 
field evolution in QuickPIC, breaks-down due to strong self- 
compression and self-steepening. 

IV. Discussion and conclusions 

Since QuickPIC simulations for such long plasmas are 
already quite large, and for some scenarios, kinetic effects 
not modeled by QuickPIC, can be relevant, we have assessed 
the possibility to use full 2D PIC simulations for quantitative 
predictions in long propagations distances. We have compared 



the results presented in Sec. Ill with full 2D simulations in 
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OSIRIS, using the same physical parameters for the laser and 
plasma channel. For the 2D OSIRIS simulation, the computa- 
tional window, that moves at the speed of light, is 230 /im wide 
and 60 /mi long. The box is divided in 400 x 2010 cells for 
the transverse and longitudinal directions, respectively, with 9 
particles per cell. 

The comparison between QuickPIC and the 2D simulations 
is shown in Fig. [5] The lower dimensionality of the 2D 
simulations lead to quantitative differences, in comparison 
with QuickPIC. In the 2D full PIC simulation, the laser 
evolution is weaker, even though the threshold power for self- 
focusing, given by P/P c = OqWq / '16\/2, is lower than in 
3D. This leads to differences of sa 30% on the laser peak 
ar> Naturally, the wakefield becomes weaker and shorter for 
the 2D OSIRIS simulations and was not even strong enough 
to self-inject plasma electrons [16]. Even so, and despite the 
quantitative differences, 2D PIC simulations and QuickPIC are 
qualitatively similar. 

These results then indicate that the numerical modeling 
of LWFA experiments, can be performed accurately only in 
3D geometries [18]. However, the use of standard full PIC 
simulations becomes prohibitive, as the computing require- 
ments increase with the interaction length associated with 
higher energy gains. Therefore, the use of QuickPIC or full 
PIC 3D simulations in boosted frames [19] (currently being 
implemented in OSIRIS [20], [21]) are valuable tools to 
examine LWFA under these conditions. 

The use of QuickPIC has allowed us to perform full-scale, 
one-to-one numerical experiments of LWFA using state-of-the- 
art experimental parameters in scenarios where full PIC 3D 
simulations are difficult to perform. We have examined the 
laser propagation in conditions which closely follow recent 
experiments [4]; for ~ 10 J laser pulses multi-GeV electron 
acceleration has been observed. 

We have shown that recent experiments have been dom- 
inated by a strong initial self-focusing. This stage is then 
followed by small spot size oscillations. In the case of the 
12 TW laser, the presence of the channel is crucial in order 
to guide the pulse after the strong initial self-focusing. In the 
25 TW laser scenario, the channel only confines the laser to 
smaller spot sizes as the laser propagates. Unlike the 12 TW 
laser, the propagation of the 25 TW laser in the pre-formed 
plasma channel, lead, after the initial self-focusing stage, to 
the continuous increase of the laser peak vector potential. 
The higher initial clq for the 25 TW laser then resulted in 
the enhancement of both longitudinal and transverse self- 
modulations, which drove a strong enough blowout region to 
guarantee guiding even in uniform plasma. 

For 10 J laser pulses, we have shown that it is possible to 
accelerate electrons up to 2.2 GeV in 1.5 cm, in an uniform 
plasma. By using a 15 J laser pulse, it is possible to generate 
at least 4 GeV electrons in 20 cm long plasma channels. The 
use of a smaller ao in this latter scenario in comparison to 
the simulation of the 10 J laser pulse, implied the use of a 
smaller background density in order to guarantee the matching 
conditions proposed by W. Lu et al [14]. Since the dephasing 
length is higher for lower densities, the maximum energy gain 
is higher in this situation. We have shown that the electron 



beam can be efficiently accelerated in long plasmas, thus 
indicating that in the near future, with the systems now coming 
online, mono-energetic electron beams can be accelerated to 
multi-GeV energies. 
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Fig. 1. (a) shows the accelerating fields at initially at s = (solid) and at 
s = 2.5 mm (dashed), after the initial self-focusing stage and (b) shows the 
laser spot size, Wo, length, to, and peak ao evolution, normalized to the initial 
values, for the 12 TW laser simulation. Both (a) and (b) refer to the situation 
where the plasma channel is present, (c) shows the comparison between the 
evolution of the spot of the laser with (solid) and without (dashed) preformed 
plasma channel. 




Fig. 2. (a) shows the accelerating fields at initially at s = (solid) and at 
s = 1.5 mm (dashed), after the initial self-focusing stage, (b) shows the laser 
spot size, Wo, length, tq, and peak ao evolution, normalized to the initial 
values for the 25 TW laser simulation. Both (a) and (b) refer to the simulations 
in the presence of the plasma channel, (c) shows the comparison between the 
evolution of the spot of the laser with (solid) and without (dashed) preformed 
plasma channel. 



: (a) s 


1 1 1 1 1 1 1 1 1 1 _ 
• — v 23 mm I 




— 12 mm - 




1 mm 


r i i i i I t -r-i- 


I |initial 



3.0 


| (b)' ' 




2.5 


% initial 23 mm 




2.0 






1.5 


% II 




1.0 


1| 1 mm / ! I 




0.5 1 


; / ! 12 mm /A J 











Fig. 3. (a) shows the evolution of the beam phase space for the 10 J laser 
simulation, where the phase rotation of the beam occurred for s = 2.3 cm. 
(b) shows the corresponding spectra revealing a quasi-mono-energetic electron 
beam with Ri 4% energy spread and with ~ 2.2 GeV after phase space 
rotation occurred. 



5 



t 1 

; (a) 


'"j"^ 2b cm" = 






^^^^ 13.3 cm ; 




6.6 cm 


f 


, 1 1 1 1 initial 



!i (b) 



3.6 cm ; 13.3 cm ■ 



20 30 



Fig. 4. (a) shows the evolution of the beam phase space for the 15 J laser 
simulation with external guiding, (b) shows the corresponding energy spectra, 
where the peak of the beam accelerated to « 4.2 GeV 
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Fig. 5. Comparisons between the 2D OSIRIS and QuickPIC simulations 
after s = 3.75 mm. (a-I) shows the on-axis longitudinal laser profile and 
(a-II) the on-axis electron density distribution for the 12TW laser pulse, (b- 
I) shows the on-axis longitudinal laser profile and (b-II) the on-axis electron 
density distribution for the 25 TW laser. 
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